Cite this Article
Cai Chen-Liang, Yu Jie, Tu Juan, Guo Xia-Sheng, Huang Pin-Tong, Zhang Dong. Interaction between encapsulated microbubbles: A finite element modelling study. Chinese Physics B, 2018, 27(8): 084302  
Permissions
Interaction between encapsulated microbubbles: A finite element modelling study
Cai Chen-Liang1, Yu Jie1, Tu Juan1, †, Guo Xia-Sheng1, Huang Pin-Tong2, Zhang Dong1, 3, ‡
Key Laboratory of Modern Acoustics (MOE), Department of Physics, Collaborative Inovation Center of Advanced Microstructure, Nanjing University, Nanjing 210093, China
Department of Ultrasound, the Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou 310009, China
The State Key Laboratory of Acoustics, Chinese Academy of Sciences, Beijing 100190, China

 

† Corresponding author. E-mail: juantu@nju.edu.cn dzhang@nju.edu.cn

Projects supported by the National Natural Science Foundation of China (Grant Nos. 11474161, 11474001, 116741731, 1774166, 11774168, 81527803, 81627802, and 81420108018), the Fundamental Research Funds for the Central Universities, China (Grant No. 020414380109, and the Qing Lan Project, China.

Abstract

Theoretical studies on the multi-bubble interaction are crucial for the in-depth understanding of the mechanism behind the applications of ultrasound contrast agents (UCAs) in clinics. A two-dimensional (2D) axisymmetric finite element model (FEM) is developed here to investigate the bubble–bubble interactions for UCAs in a fluidic environment. The effect of the driving frequency and the bubble size on the bubble interaction tendency (viz., bubbles’ attraction and repulsion), as well as the influences of bubble shell mechanical parameters (viz., surface tension coefficient and viscosity coefficient) are discussed. Based on FEM simulations, the temporal evolution of the bubbles’ radii, the bubble–bubble distance, and the distribution of the velocity field in the surrounding fluid are investigated in detail. The results suggest that for the interacting bubble–bubble couple, the overall translational tendency should be determined by the relationship between the driving frequency and their resonance frequencies. When the driving frequency falls between the resonance frequencies of two bubbles with different sizes, they will repel each other, otherwise they will attract each other. For constant acoustic driving parameters used in this paper, the changing rate of the bubble radius decreases as the viscosity coefficient increases, and increases first then decreases as the bubble shell surface tension coefficient increases, which means that the strength of bubble–bubble interaction could be adjusted by changing the bubble shell visco-elasticity coefficients. The current work should provide a powerful explanation for the accumulation observations in an experiment, and provide a fundamental theoretical support for the applications of UCAs in clinics.

PACS: 43.25.Yw;43.35.Wa;43.80.+p
Keyword:ultrasound contrast agent microbubbles;bubble-bubble interaction;finite element model;shell parameter
1. Introduction

Ultrasound contrast agents (UCAs) are encapsulated microbubbles filled with gas, which have been initially utilized in clinics as contrast agents for ultrasound diagnostic imaging because of their higher echogenicity than those of the background tissues.[1,2] Recently, the UCAs have shown growing potential in therapeutic ultrasound applications, in which ultrasound bio-effects can be enhanced due to the cavitation of microbubbles, arousing the increasing interest for targeted drug delivery, high intensity focused ultrasound (HIFU), and other applications.[3–10]

For better understanding the mechanisms behind the UCA applications in diagnostic and therapeutic applications, an increasing number of researches have been conducted on the dynamics of encapsulated microbubbles. The existing studies modified the classic vibration equation of encapsulated bubbles, which is known as the Rayleigh–Plesset equation,[11] by considering the effects of shell properties, such as surface tension and dilatational viscosity. Moreover, the interaction among several encapsulated bubbles, and the kinetic behavior of UCAs inside biological tissues were also investigated to systematically simulate the dynamics of microbubbles in vivo. Doinikov et al. conducted a preliminary research on the bubble–bubble interaction based on the numerical algorithm and indicated that when exposed to strong acoustic pressure, microbubbles can form a bound pair with a steady spacing rather than collide and coalesce.[12] Martynov et al. studied numerically the effect of the elastic wall of vessels on the oscillation of microbubbles based on the finite element method (FEM), finding that the bubble oscillations in a finite-length vessel were characterized by a spectrum of frequencies, with distinguishable high-frequency and low-frequency modes.[13] Qin et al. established a numerical model of the interaction between the bubble and the vessel wall based on the FEM, demonstrating that the vibration of the bubble could increase the permeability of blood vessels.[4,14] In an earlier research, Miao and Gracewski developed a two-dimensional (2D) hybrid model combining the FEM and the boundary element method (BEM) algorithm for studying the interaction between a microbubble and either a deformable structure or a micro-vessel. They found that the presence of a tube mitigated the bubble expansion.[15] In addition, the maximum tube dilation and the maximum tensile hoop stress occurred during the tensile portion of the acoustic excitation, well before the maximum bubble. Chen et al. performed numerical simulations based on a 2D asymmetric finite element model to investigate the influence of both acoustic driving parameters and material properties on the dynamic interaction in the bubble–blood–vessel system, and the results indicated that the constrained effect of the blood vessel along the radial direction would induce the asymmetric bubble oscillation and vessel deformation, as well as shifting the bubble resonance frequency toward the higher frequency range.[16]

Although the forces exerted on microbubbles induced by bubble–bubble interactions can be calculated numerically, there is a lack of detailed information about the velocity field in surrounding fluid and bubbles’ relative motions during the interaction. In this paper, an FEM model regarding the dynamic behaviors of bubble–bubble interaction is developed, in which the modification of the bubble dynamic model proposed by Chatterjee and Sarkar is considered.[17] Based on the FEM model, the interactions between bubbles are systematically investigated, including the changes in the bubble radius and relative bubble–bubble distance during the bubble–bubble interaction, as well as the distribution of the velocity field in the surrounding fluid, which plays a significant role in microbubble dynamics. In addition, the effect of bubble shell mechanical parameters (viz., surface tension coefficient and viscosity coefficient) on bubble–bubble interaction are also investigated.

2. Hydrodynamics of bubble dynamics

The dynamics of microbubbles in an infinite liquid environment is described by the coupling motion of the gas-liquid boundary. When a microbubble is exposed to ultrasound, the fluid around the bubble satisfies two basic equations, i.e., the continuity equation and the Navier-Stokes equation:

where u denotes the velocity, p is the pressure, I is the unit tensor matrix, ρl is the fluid density, and μl is the dynamic viscosity. The assumption of incompressible fluid leads to the fact that the liquid should meet the state equation:

In addition, the gas inside the bubble should satisfy the state equation of the ideal gas, since the bubble vibration is considered to be an adiabatic process:

where pg is the gas pressure inside the bubble with an initial value of pg0, V is the volume of the bubble, V0 is the initial volume, and γ is the polytropic exponent.

For the sake of simplifying the numerical simulation, only the dynamic responses induced by the bubble–bubble interactions are considered here. Figure 1 illustrates the bubble–bubble interaction diagram, where D(t) is the instantaneous distance between the centers of the two bubbles, and Rj(t) (j = 1 or 2) is the instantaneous radius of an individual bubble. In order to maintain the spherical vibration of the bubble, Rj is set to be smaller enough than D(t), and the ratio of D(t) to Rj(t) is larger than 10 in the simulation.[18,19]

Fig. 1. (color online) Geometric illustration of bubble–bubble interaction.

In addition, the common symmetric axis of two bubbles is defined as the z axis in this paper and the initial midpoint of the line segment between two bubbles’ centers is set to be as the original point (viz., point O).

If the bubbles are considered as air-free bubbles, the gas–liquid boundary should satisfy the following boundary conditions:

Given that the bubbles are encapsulated bubbles, the equation will be modified by taking into account the influence of the shell properties on the bubble dynamics as follows:[17]
where σg is the surface tension of the air bubble, while σe is the surface tension coefficient of the encapsulated bubble, and κs is the viscosity coefficient of the shell.

In addition, the pressure at infinity is required to satisfy the boundary condition:

where p0 is the hydrostatic pressure of the environment, and pi(t) is the driving pressure of ultrasound.

Considering that the gas–liquid boundary is initially balanced, the initial condition needs to be satisfied[20]

Based on the above equations, the linear resonance frequency of air bubble can be calculated from the following equation:

3. Finite element simulation

Based on the computational fluid dynamics (CFD) module and the Arbitrary Lagrangian–Eulerian (ALE) method, an FEM model for simulating the interactions between two encapsulated bubbles is established with Comsol Multiphysics 5.2a. To simplify the computation, the system is treated as being 2D axisymmetric. Figure 2 shows the schematic of the model and the mesh configuration at the bubble boundary. The fluid domain is decomposed into triangular meshes, and the mesh density of the bubble–fluid interface is enhanced to assure that the deformation of the bubble shell can be resolved more precisely. On the spherical bubble–fluid boundary, the mesh element size is set to be smaller than 1/10 of the bubble radius, and the rest of the fluid domain was decomposed with the self-adaption method of the FEM software. Some physical parameters used in the simulation are listed as follows:[21] density of fluid ρl = 103 kg/m3, dynamic viscosity of fluid μl = 10−3 Pa·s, gas polytropic exponent γ = 1.07, surface tension of shell σe = 0.32 N/m, viscosity coefficient of shell κs = 4 × 10−9 kg/s, hydrostatic pressure of environment p0 = 1.013 × 105 Pa, and initial bubble–bubble distance D(0) = 30 μm, which means that bubble 1 and bubble 2 are initially located at z10 = −15 μm and z20 = 15 μm, respectively.

Fig. 2. Schematic and meshes of the model in Comsol.
4. Results and discussion

In the simulation, the driving ultrasound pressure is set to be pi(t) = pmax · sin(2πft), where f denotes the driving frequency and pmax refers to the amplitude of pressure, which is typically set to be pmax = 0.05 MPa.

4.1. Translational phenomenon of bubble–bubble interaction
4.1.1. Interaction between microbubbles with the same size

In this part, the initial radii of two bubbles are set to be R1(0) = R2(0) = 1.5 μm, and the resonant frequencies of the two bubbles can be calculated by using Eq. (9) to be f10 = f20 = 3.7 MHz. Thus, the driving frequency of ultrasound is set to be the same as the linear resonance frequency of the bubbles, in order to observe the bubble interaction more clearly. The acoustic driving pressure is pmax = 0.05 MPa.

Figure 3 shows the temporal variation of the bubble radius and the bubble–bubble distance. It is indicated that the bubbles interact with each other through the coupling of the fluid, resulting in the change of bubble–bubble distance, which decreases from 30 μm to 29.81 μm. The results indicate that the overall trend of the translational motion induced by the interaction between two equal-sized bubbles is a process of mutual attraction. In addition, an obvious oscillation can be observed in each excitation cycle, indicating that attractive force and repulsive force alternate between the bubbles during the interaction, which might be closely related to the dynamic motions of the velocity field in the surrounding fluid. The attraction can be explained by the acoustic radiation force, which consists of two components, i.e., a primary force (directing away from the source), and a secondary force (typically attractive between bubbles).[22] As indicated by Rychak et al., the secondary radiation force between two identical microbubbles equidistant from the acoustic source could cause the bubbles to attract each other, forming the specific target accumulation of microbubbles.[23]

Fig. 3. Temporal evolution of (a) the radius and (b) the distance between bubbles with the same size. The part enclosed by a dashed circle is enlarged as shown in the inset. The initial bubble radii are R1(0) = R2(0) = 1.5 μm. The acoustic driving parameters are f = 3.7 MHz and pmax = 0.05 MPa.
4.1.2. Interaction between encapsulated microbubbles with different sizes

In clinical practice, it is difficult to produce microbubbles with an identical size precisely. Thus, more detailed studies are performed to simulate the interactions between two bubbles with different sizes. Here the initial radii of the bubbles are set to be R1(0) = 1.5 μm and R2(0) = 4 μm, respectively, and the initial bubble–bubble distance is set to be D(0) = 30 μm. Consequently, the resonance frequencies of the two bubbles are calculated to be f10 = 3.7 MHz and f20 = 1.0 MHz.

Figures 4–6 show the curves of the bubble radii (a) and the bubble–bubble distance (b) versus time when the driving frequency is f = 2.4 MHz (i.e., f10 < f < f20), 4.7 MHz (i.e., f > max (f10, f20)) and 0.7 MHz (i.e., f < min(f10, f20)), respectively, and the acoustic pressure amplitude is kept to be pmax = 0.05 MPa. As shown in Fig. 4, when the driving frequency is between the resonant frequencies of two bubbles (i.e., f = 2.4 MHz), generally a repulsive motion of bubbles is observed for interacted bubbles. When the driving frequency is greater or smaller than the resonant frequency of the bubbles (i.e., f = 4.7 MHz or f = 0.7 MHz), the overall interaction motion demonstrated by these two bubbles is to attract each other as shown in Figs. 5 and 6, which is consistent with measured results of the bubble translation in a multi-bubble environment conducted by Xi.[24]

Fig. 4. Temporal evolution of the two bubble radii (see (a) and (b), respectively), and (c) bubble–bubble distance. Inset shows a magnified part of the curve, enveloped in a dashed circle. Initial bubble radii are R1(0) = 1.5 μm and R2(0) = 4.0 μm, respectively. Acoustic driving parameters are f = 2.4 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.
Fig. 5. Temporal evolution of the two bubble radii (see (a) and (b), respectively), and (c) bubble–bubble distance. Inset shows part of the curve, enveloped in a dashed circle. Initial bubble radii are R1(0) = 1.5 μm and R2(0) = 4.0 μm, respectively. Acoustic driving parameters are f = 4.7 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.
Fig. 6. Temporal evolution of the two bubble radii (see (a) and (b), respectively), and (c) bubble–bubble distance. Inset shows a magnified part of the curve, enveloped in a dashed circle. The initial bubble radii are R1(0) = 1.5 μm and R2(0) = 4.0 μm, respectively. The acoustic driving parameters are f = 0.7 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.

Although there is an overall translational trend of repulsion or attraction between two bubbles, the bubble–bubble distances always exhibit obvious oscillating motions in these three cases (Figs. 4–6), which is similar to the case of two bubbles with equal size (Fig. 3(b)). This phenomenon should be attributed to the gradient variation in the velocity field of the fluid surrounding the bubbles, and will be discussed in the following subsection.

It is worth mentioning that the two bubbles with different sizes exhibit different behaviors when they are exposed to ultrasound at a fixed frequency. For example, bubble 1 with a smaller size exhibits a larger oscillation amplitude at f = 4.7 MHz as its resonance frequency is closer to the driving frequency (Figs. 5(a) and 5(b)), whereas bubble 2 with a larger size shows a more significant change in the radius when driven by ultrasound at a frequency of f = 0.7 MHz (Figs. 6(a) and 6(b)).

4.2. Velocity field of fluid around bubble

As illustrated in Figs. 3–6, there is an identical feature in the four D(t) curves under four different conditions that the instantaneous positions of bubbles’ centers show quasi-periodic vibration during bubble–bubble interaction, no matter whether the overall trend is repulsive or attractive, which is difficult to explain with the classical Bjerknes force. Here, the instantaneous distributions of the velocity in the fluid surrounding the bubbles are simulated to explore the possible mechanism behind this phenomenon.

In order to obtain an in-depth understanding of the influence of the fluid velocity field on the motion of the bubbles, two moments during bubble–bubble interactions under ultrasound exposures are selected for individual cases, when two bubbles are instantaneously repulsive and attractive, respectively. Figure 7 demonstrates the velocity field in the fluid around the bubble–bubble couple with the same size (e.g., R1(0) = R2(0) = 1.5 μm) at two instantaneous time points, viz., t1 = 15.38 μs (Fig. 7(a), instantaneously repelling motion) and t2 = 16.60 μs (Fig. 7(b), instantaneously attracting motion). The details of the instantaneous overall translational trend and a single bubble motion at these two time points are illustrated in the inset of Fig. 3(b) and Fig. 7(c), respectively. It is obvious that the velocity field around two bubbles with the same size is always symmetric and the dynamic motions of these two bubbles are exactly the same, so that the two bubbles contribute equally during their interacting process.

Fig. 7. (color online) Distribution of velocity field in fluid at (a) t1 = 15.38 μs and (b) t2 = 16.60 μs around two bubbles with the same size, and (c) instantaneous |z|–t curve. Initial bubble radii are R1(0) = R2(0) = 1.5 μm. The acoustic driving parameters are f = 3.7 MHz and pmax = 0.05 MPa.

However, for the bubbles with different radii (e.g., R1(0) = 1.5 μm and R2(0) = 4.0 μm), an asymmetric velocity field can be observed around these two bubbles and different oscillating motions are demonstrated by individual bubbles. For instance, when the driving frequency of ultrasound is 2.4 MHz (i.e., f10 < f < f20), the overall translational trend between the two bubbles with different sizes is repulsive, although the whole repelling process is a quasi-periodic oscillation process, as shown in Fig. 4(c). Figure 8 shows the velocity field in the fluid around the bubbles at two instantaneous time points (viz., t1 = 18.35 μs and t2 = 18.55 μs). The corresponding instantaneous interacting motions at these two points are repulsive and attractive respectively as illustrated in the inset in Fig. 4(c). It can be noticed in Fig. 8(a), at t1 = 18.35 μs, the larger bubble 2 is in the state of expansion and the smaller bubble 1 is in the state of contraction, so that the velocity field between these two bubbles points from bubble 2 to bubble 1. In terms of the smaller bubble 1, the inside velocity field is stronger than the outside velocity field due to the superposition with the velocity field generated by the expansion of the larger bubble 2, so that bubble 1 is pushed away from bubble 2. In terms of bubble 2, due to the superposition with the velocity field generated by the contraction of bubble 1, the inside velocity field overweighs the outside velocity field and bubble 2 moves towards bubble 1. In general, through the fluid as a medium, bubble 2 obtains an attractive force while bubble 1 obtains a repulsive force, because they oscillate out of phase. At t2 = 18.55 μs, the larger bubble 2 is contracting and the smaller bubble 1 is expanding, so the direction of the velocity field in the fluid flow points from the bubble 1 to the bubble 2, resulting in the enhancement of the velocity field between the two bubbles, and thus the bubble 2 is pushed away from the original point and the bubble 1 is pulled towards the original point.

Fig. 8. (color online) Distribution of velocity field in the fluid at (a) t1 = 18.35 μs, (b) t2 = 18.55 μs, and (c) instantaneous |z|–t curve. Acoustic driving parameters are f = 2.4 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.

When the driving ultrasound frequency is set to be f = 4.7 MHz (i.e., f > max (f10, f20)), the interaction between two bubbles shows an overall tendency of attraction with oscillation. Figure 9 shows the fluid velocity field distributions around the two bubbles at two different time points, t1 = 16.15 μs and t2 = 17.55 μs. Since these two bubbles always expand or contract simultaneously, their velocity fields counteract each other, which makes their instantaneous moving motions dominated by their outside velocity fields. Therefore, it is observed that at t1 = 16.15 μs, the direction of the outside velocity field in the fluid points from two bubbles, which makes the bubble–bubble couple mutually repulsive. On the contrary, at t2 = 17.55 μs, the outside velocity field points into both bubbles, which makes them attract each other at this moment. This tendency of attraction and repulsion can also be observed from the curve D(t) at the same time (see the inset in Fig. 5(c)).

Fig. 9. (color online) Distribution of velocity field in fluid at (a) t1 = 16.15 μs and (b) t2 = 17.55 μs, and (c) instantaneous |z|–t curve. The acoustic driving parameters are f = 4.7 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.

When the driving frequency of ultrasound is set to be f = 0.7 MHz (i.e., f < min (f10,f20), the interaction between two bubbles shows an overall tendency of mutual attraction with oscillation. Figure 10 shows the velocity field distributions around two bubbles at two different time points, t1 = 15.85 μs and t2 = 16.30 μs. As shown in Fig. 10(a), at t1 = 15.85 μs, it is noticed that these two bubbles are both in the state of contraction. Therefore, their velocity fields counteract each other, causing their instantaneous translational motions to be dominated by their outside velocity fields and be pushed closer towards each other. As shown in Fig. 10(b), at t2 = 16.30 μs, the velocity generated by the expansion of the larger bubble 2 exerts a repulsive force on bubble 1, causing bubble 1 to move away from bubble 2. Meanwhile, due to the counteraction of the velocity generated by bubble 1, the velocity field inside bubble 2 is also weaker than its outside velocity field, so that it moves away from bubble 1. This tendency of attraction and repulsion could be confirmed from the curve D(t) at the same time points (see the inset in Fig. 6(c)).

Fig. 10. (color online) Distribution of velocity field in fluid at (a) t1 = 15.85 μs, (b) t2 = 16.30 μs, and (c) instantaneous |z|–t curve. The acoustic driving parameters are f = 0.7 MHz, pmax = 0.05 MPa, f10 = 3.7 MHz, and f20 = 1.0 MHz.

The above observations suggest that for two bubbles with different sizes, the overall translational trend (i.e., repulsion or attraction) should be determined by the relationship between the driving frequency and the resonance frequencies of the two bubbles, which is consistent with the numerical investigation done by Pelekasis and Tsamopoulus in the motion of two bubbles induced by an oscillatory disturbance at the ambient pressure.[25] Generally speaking, when the driving frequency is between the resonant frequencies of these two bubbles (e.g., f = 2.4 MHz), the bubble–bubble couple demonstrates a repulsive motion (Fig. 4(c)). However, if the driving frequency is greater or smaller than the resonant frequency of both bubbles (e.g., f = 4.7 MHz or f = 0.7 MHz), the bubble–bubble couple will exhibit an attraction tendency (Figs. 5(c) and 6(c)).

In addition, the details of the fluid velocity field around the bubbles show that there are always periodic oscillations in the variation of bubble–bubble distance, and the instantaneous motions of these two bubbles depend on the surrounding velocity field. If the velocity field superposes over each other between these two bubbles, the inside velocity field overweighs the outside velocity field and the instantaneous translational motion is dominated by the inside velocity field, and vice versa. Fundamentally, it is the motion phase that determines the superposition and counteraction of the velocity fields between these two bubbles. Here, with the FEM model, it is convenient to derive some segments or points in the time sequence to observe the detail of the instantaneous flow velocity field and vibration motion of the bubbles. Thus, one can see from Fig. 8 that when two bubbles with different sizes are exposed to ultrasound whose frequency falls in the interval between their resonance frequencies, their oscillating motion is always out of phase, almost having a difference of π, which makes their velocity fields superposed and the instantaneous translational motion is dominated by the inside velocity field. Otherwise, when the driven frequency is greater or smaller than their resonance frequency, the two bubbles always expand or contract simultaneously to exert a counteracting effect on the velocity fields between these two bubbles, so that their instantaneous translational motion is dominated by their outside velocity field as shown in Figs. 9 and 10. It can help to explain the phenomenon observed by Masato that the bubbles generally vibrating in-phase are easy to attract each other, while repulsive motion is usually observed for those nearly vibrating out of phase.[26]

4.3. Influence of shell parameters on bubble interaction

It has also been well accepted that the presence of the microbubble shell material would play an influencing role in the microbubble dynamic response.[27–37] Therefore, the influence of UCA shell visco-elasticity properties (viz., surface tension coefficient and viscosity coefficient) on the bubble–bubble interaction is also studied in the following.

For simplicity, only the interactions between bubbles with the same size are considered in the present work. The radii of the two bubbles are set to be R1(0) = R2(0) = 1.5 μm, and the ultrasound driving frequency is f = 3.7 MHz, with a pressure amplitude of pmax = 0.05 MPa. Based on the multiphysics finite element modeling with Comsol 5.2a, one can obtain the maximum bubble radius Rmax, the minimum bubble radius Rmin, and the value of D(t) at any typical time point (e.g., t = 100 μs).

Figure 11 shows the changes in the bubble radii and bubble–bubble distance as the bubble shell viscosity coefficient increases, where the bubble shell surface tension coefficient is fixed to be 0.32 N/m. It is found that the changing rate of the bubble radius becomes smaller as the viscosity coefficient increases, which indicates that the bubble response to the ultrasound is weakened due to the viscous damping. Consequently, the changing rate of the velocity gradient of the fluid around the bubble is reduced, resulting in the decrease of the change of the bubble–bubble distance of the two bubbles. As discussed by Sadighi-Bonabi et al., in clinical applications, the secondary Bjerknes force could be reduced by increasing the fluid viscosity in the surrounding medium.[38]

Fig. 11. Dependence of (a) Rmax, (b) Rmin, and (c) D(t) on the bubble shell viscosity coefficient κs, with R1(0) = R2(0) = 1.5 μm, pmax = 0.05 MPa, f = 3.7 MHz, t = 100 μs, and bubble shell surface tension coefficient being fixed to be 0.32 N/m.

Then, with the bubble shell viscosity coefficient kept at a constant value κs = 4 × 10−9 kg/s, the influence of bubble shell surface tension coefficient is investigated. Figure 12 shows the dependence of the bubble radii and bubble–bubble distance on the bubble shell surface tension coefficient. It is observed that all the variables achieve the peaks at σe = 0.32 N/m. The reason for this phenomenon is that the resonance frequency of the bubble is close to the driving frequency when the bubble shell surface tension coefficient is close to 0.32 N/m, so that the response of bubble to the external ultrasound could be maximized due to the resonance effect, which results in the maximum velocity gradient in the fluid around the bubbles, and the change of the bubble–bubble distance is also maximized. Based on the resonance phenomenon, it comes to a conclusion that when the acoustic driving parameters are fixed in an experiment, the strength of bubble–bubble interaction could be adjusted by changing the bubble shell visco-elasticity coefficient.

Fig. 12. Dependence of (a) Rmax, (b) Rmin, and (c) D(t) on the bubble shell surface tension coefficient σe, with R1(0) = R2(0) = 1.5 μm, pmax = 0.05 MPa, f = 3.7 MHz, t = 100 μs, and bubble shell viscosity coefficient kept at a constant value κs = 4 × 10−9 kg/s.
5. Conclusions

In this paper, a model of microbubble vibration and bubble–bubble interaction is built from the perspective of fluid mechanics based on the multi-physics finite element method. Compared with the traditional numerical method, the FEM has the advantage to obtain a 2D instantaneous velocity field distribution diagram, which provides the possibility to explain the interaction mechanism of bubble interaction from the viewpoint of fluid mechanics, which acts as a foundation for the in-depth study of the microbubble dynamics by studying the relationship between the bubble vibrations, the interaction between bubbles, and the visco-elastic properties of the bubble shell. The results indicate that for the interacting bubbles, the overall translational trend is dominated by the relationship between the driving frequency and their resonance frequency. For constant acoustic driving parameters, the bubble–bubble interactions can be regulated by adjusting the shell visco-elasticity coefficients of the bubbles. Typically, the interaction strength decreases as the viscosity coefficient increases, and increases first then decreases as the bubble shell surface tension coefficient increases.

The definition of the microbubble shell parameters in the finite element model in this paper is based on the modified microbubble model proposed by Sarkar.[19] This model neglects the variation of the film thickness and the changing of the shell parameter with the radius of bubble. In the following work, the modified parameters of Marmottant, Doinkov, and other more complex models can be modified to further improve the accuracy of the FEM simulation results. In this model, we assume that the bubble is in an infinite space, which may be further extended to more complex conditions. In future work, we can consider the interaction of multiple bubbles and the constraint of the wall of the blood vessel on the bubble vibration, and thus making the simulation closer to the actual experimental situation. In addition, considering the convergence of the finite element calculation and the mesh distortion, the intensity of ultrasonic excitation is relatively small, which is the limitation of the finite element method. To consider the interaction between the bubbles driven by ultrasound of relatively large pressure amplitude, the finite element method could be combined with the boundary element method, which is a subject that can be explored in the future.

Reference
[1] Zhang D Gong Y Gong X Liu Z Tan K Zheng H 2007 Phys. Med. Biol. 52 5531
[2] Tu J Guan J F Matula T J Crum L A Wei R J 2008 Chin. Phys. Lett. 25 172
[3] Ferrara K Pollard R Borden M 2007 Ann. Rev. Biomed. Eng. 9 415
[4] Qin S Caskey C F Ferrara K W 2009 Phys. Med. Biol. 54 R27
[5] Stride E P Coussios C C 2010 Proc. Inst. Mech. Eng. Part. H J. Eng. Med. 224 171
[6] Zhou Y F 2011 World J. Clin. Oncol. 2 8
[7] Cai X W Yang F Gu 2012 Theranostics 2 103
[8] Arvanitis C D Gregory T C Nathan M 2015 IEEE Trans. Med. Imaging 34 1270
[9] Lee J Min H S You D G Kim K Kwon I C Rhim T Lee K Y 2016 J. Control. Release 223 197
[10] Zhou Y F Gao X W 2016 Phys. Med. Biol. 61 6651
[11] Plesset M S Prosperetti A 1977 Annu. Rev. Fluid Mech. 9 145
[12] Doinikov A A 2001 Phys. Rev. 64 026301
[13] Martynov S Stride E Saffari N 2009 J. Acoust. Soc. Am. 126 2963
[14] Qin S Ferrara K W 2006 Phys. Med. Biol. 51 5065
[15] Miao H Y Gracewski S M 2008 Comput. Mech. 42 95
[16] Chen C Gu Y Tu J Guo X Zhang D 2016 Ultrasonics 66 54
[17] Chatterjee D Sarkar K 2003 Ultrasound Med. Biol. 29 1749
[18] Fuster D Conoir J M Colonius T 2014 Phys. Rev. 90 063010
[19] Cui B Ni B Wu Q 2016 Adv. Mech. Eng. 8 1687814016631708
[20] Leighton T 1994 The Acoustic Bubble London Academic Press 10.1016/B978-0-12-441920-9.X5001-9
[21] Tu J Guan J Qiu Y Matula T J 2009 J. Acous. Soc. Am. 126 2954
[22] Dayton P A Morgan K E Klibanov A L Brandenburger G Nightingale K R Ferrara K W 1997 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 44 1264
[23] Rychak J J Klibanov A L Hossack J A 2005 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 421
[24] Xi X 2013 Controlled Translation and Oscillation of Micro-bubbles Near a Surface in an Acoustic Standing Wave Field Ph. D. Diserlation London Imperial College
[25] Pelekasis N A Tsamopoulos J 1993 J. Fluid Mech. 254 501
[26] Masato I 2003 Phys. Rev. 67 056617
[27] de Jong N Hoff L 1993 Ultrasonics 31 175
[28] Shi W T Forsberg F Hall A L Chiao R Y Liu J B Miller S Thomenius K E Wheatley M A Goldberg B B 1999 Ultrasonic Imaging 21 79
[29] Frinking P J de Jong N 1998 Ultrasound Med. Biol. 24 523
[30] Guo X Li Q Zhang Z Zhang D Tu J 2013 J. Acoust. Soc. Am. 134 1622
[31] Guo G Lu L Yin L Tu J Guo X Wu J Xu D Zhang D 2014 Phys. Med. Biol. 59 6729
[32] Wang L Tu J Guo X Xu D Zhang D 2014 Chin. Phys. 23 124302
[33] Goldberg B B Raichlen J S Forsberg F 2001 Ultrasound contrast agents: basic principles and clinical applications Boca Raton CRC Press
[34] Gorce J M Arditi M Schneider M 2000 Invest. Radiol. 35 661
[35] Tu J Swalwell J E Giraud D Cui W Chen W Matula T J 2011 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58 955
[36] Doinikov A A Bouakaz A 2011 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58 981
[37] Doinikov A A Bouakaz A 2013 Phys. Med. Biol. 58 6797
[38] Sadighi-Bonabi R Rezaee N Ebrahimi H Mirheydari M 2010 Phys. Rev. 82 016316